#include "Collision1D_Ionization.h"
#include "../Parallel/ParallelMPI.h"
#include "../Field/Field2D.h"


#include <cmath>
#include <iomanip>
#include <algorithm>
#include <ostream>
#include <sstream>

using namespace std;

Collision1D_Ionization::Collision1D_Ionization(PicParams* params_, CollisionParameter collision_param_, int collision_number_) :
    Collision1D(params_, collision_param_, collision_number_)
{
    //cout<<"collision111"<<endl;
    read_cross_section();
    //cout<<"collision222"<<endl;
    energy_threshold = crossSection[0][0];
    //cout<<"collision333"<<endl;
}

Collision1D_Ionization::~Collision1D_Ionization()
{

}


void Collision1D_Ionization::collide(PhysicalField* fields, vector<Species*>& vecSpecies, Diagnostic* diag, int itime)
{
    vector<int> index1, index2;
    vector<int> n1, n2;
    int n1_max, n2_max;
    vector<double> density1, density2;
    double density1_max, density2_max;
    vector<double> velocity_unit(3, 0.0), velocity_temp(3, 0.0);
    int collide_times = 0;
    int npairs_max;
    double npairs_double_max;
    unsigned int i1, i2, i3;
    double m1, m2, m3;

    double  sigma_cr, sigma_cr_max, niu, niu_max, ke1, ke_primary, ke_secondary,
            ran, P_collision;
    double  v_square, v_magnitude, v_magnitude_primary, v_magnitude_secondary;

    int i_space_global;
    double ke_radiative;


    if(collision_param.timestep_collision == 0 || itime % collision_param.timestep_collision != 0)
    {
        return;
    }


    Diagnostic1D *diag1D = static_cast<Diagnostic1D*>(diag);

    //electons
    auto s1 = static_cast<Species1D*>(vecSpecies[collision_param.species_number_list[0]]);

    //atom
    auto s2 = static_cast<Species1D*>(vecSpecies[collision_param.species_number_list[1]]);

    //ion
    auto s3 = static_cast<Species1D*>(vecSpecies[collision_param.species_number_list[2]]);

    m1 = s1->species_param.mass;
    m2 = s2->species_param.mass;
    m3 = s3->species_param.mass;

    int n_space = params->dims_space[0];

    n1.resize(n_space);
    n2.resize(n_space);
    density2.resize(n_space);

    density2_max = 0.0;
    n1_max = 0;

    for(int i_space = 0; i_space < n_space; i_space++)
    {
        n1[i_space] = s1->particles_[i_space]->size();
        n2[i_space] = s2->particles_[i_space]->size();

        density2[i_space] = n2[i_space] * s2->species_param.weight;

        if(density2[i_space] > density2_max) { density2_max = density2[i_space]; }
        if(n1[i_space] > n1_max) { n1_max = n1[i_space]; }
    }

    collide_times = 0;
    sigma_cr_max = max_cv(s1);
    //niu_max = sigma_cr_max * density2_max;

    //double npairs_double_max_coef = (1.0 - exp(-niu_max * collision_param.timestep_collision * params->dt * collision_param.time_zoom_factor));

    for(int i_space = 0; i_space < n_space ; i_space++)
    {
        //npairs_double_max = n1[i_space] * (1.0 - exp(-niu_max * collision_param.timestep_collision * params->dt * collision_param.time_zoom_factor));
        npairs_double_max = n1[i_space] * (1.0 - exp(-density2[i_space] * sigma_cr_max * collision_param.timestep_collision * params->dt));
        npairs_double_max *=  collision_param.time_zoom_factor;

        //determine if zoom collision frequency
        if(collision_param.region_collision_zoom_start.size() > 0 && collision_param.region_collision_zoom_end.size() > 0)
        {
            if(params->local_min[0] + (i_space + 1) * params->cell_length[0] < collision_param.region_collision_zoom_start[0]
               || params->local_min[0] + i_space * params->cell_length[0] > collision_param.region_collision_zoom_end[0])
            {
                npairs_double_max /= collision_param.time_zoom_factor;
            }
        }


        if(n1[i_space] == 0 || n2[i_space] == 0)
        {
            continue;
        }
        
        auto cell_particles1 = s1->particles_[i_space];
        auto cell_particles2 = s2->particles_[i_space];
        auto cell_particles3 = s3->particles_[i_space];

        index1.resize(n1[i_space]);
        for(int i = 0; i < index1.size(); i++)
        {
            index1[i] = i;
        }
        random_shuffle(index1.begin(), index1.end());

        index2.resize(n2[i_space]);
        for(int i = 0; i < index2.size(); i++)
        {
            index2[i] = i;
        }
        random_shuffle(index2.begin(), index2.end());

        npairs_double_[i_space] = npairs_double_max;

        npairs_[i_space] = npairs_double_[i_space];

        //See equations in http://dx.doi.org/10.1063/1.4742167
        npairs_rem_[i_space] += ( npairs_double_[i_space] - npairs_[i_space] );
        if(npairs_rem_[i_space] >= 1.0)
        {
            npairs_rem_[i_space] = npairs_rem_[i_space] - 1.0;
            npairs_[i_space]++;
        }

        if(npairs_[i_space] > n1[i_space])
        {
            cout<<"npairs is larger than the particle number in a cell!!!"<<endl;
            cout<<"npairs, n1 are: "<<npairs_[i_space]<<" "<<n1[i_space]<<endl;
            npairs_[i_space] = n1[i_space];
        }
        if(npairs_[i_space] > n2[i_space])
        {
            cout<<"npairs is larger than the particle number in a cell!!!"<<endl;
            cout<<"npairs, n2 are: "<<npairs_[i_space]<<" "<<n2[i_space]<<endl;
            npairs_[i_space] = n2[i_space];
        }


        for(int i = 0; i < npairs_[i_space]; i++)
        {
            if(npairs_[i_space] > 1)
            {
                cout<<"npairs_[i_space] = "<<npairs_[i_space]<<endl;
            }
            //MESSAGE("nparis111"<<"  "<<i);
            i1 = index1[i];
            i2 = index2[i];
            auto& p1 = cell_particles1->data_[i1];
            auto& p2 = cell_particles2->data_[i2];

            v_square = pow(p1.velocity[0],2) + pow(p1.velocity[1],2) + pow(p1.velocity[2],2);
            v_magnitude = sqrt(v_square);
            //kinetic energy of species1 (incident electrons)
            ke1 = 0.5 * m1 * v_square;
            ke_primary = ke1 - energy_threshold * params->const_e;

            //the energy of the secondary electron
            ran = (double)rand() / RAND_MAX;
            ke_secondary = 10.0 * tan(ran * atan( (ke_primary/params->const_e) / 20.0));
            ke_secondary *= params->const_e;
            //ke_secondary = 0.5 * ke_primary;
            //the energy of the primary electron
            ke_primary -= ke_secondary;
            v_magnitude_primary = sqrt( 2.0 * ke_primary / m1 );
            v_magnitude_secondary = sqrt( 2.0 * ke_secondary / m1 );

            sigma_cr = v_magnitude * interplate_cross_section( ke1 / params->const_e );
            //niu = density2[i_space] * sigma_cr;
            //P_collision = niu / niu_max;
            P_collision = sigma_cr / sigma_cr_max;
            //Generate a random number between 0 and 1
            double ran_p = (double)rand() / RAND_MAX;

            if(ran_p < P_collision)
            {
                //cout<<"collide happen "<<endl;
                cell_particles2->add_particle_delete_list(i2);

                //calculate the scatter velocity of primary electron
                velocity_unit[0] = p1.velocity[0] / v_magnitude;
                velocity_unit[1] = p1.velocity[1] / v_magnitude;
                velocity_unit[2] = p1.velocity[2] / v_magnitude;
                //MESSAGE("v_magnitude"<<"  "<<v_magnitude_primary<<"  "<<v_magnitude_secondary);
                //MESSAGE("velocity1"<<" "<<p1->velocity(0, i1)<<"  "<<p1->velocity(1, i1)<<"  "<<p1->velocity(2, i1));
                calculate_scatter_velocity(ke_primary/params->const_e, v_magnitude_primary, m1, m2, velocity_unit, velocity_temp);
                p1.velocity[0] = velocity_temp[0];
                p1.velocity[1] = velocity_temp[1];
                p1.velocity[2] = velocity_temp[2];

                //calculate the scatter velocity of secondary electron
                calculate_scatter_velocity(ke_secondary/params->const_e, v_magnitude_secondary, m1, m2, velocity_unit, velocity_temp);
                //create new particle in the end of p1, we should sort_part when all bins are done!!!
                cell_particles1->add_particles(1);

                int id_new = cell_particles1->size() - 1;
                auto& p1_new = cell_particles1->data_[id_new];
                p1_new.velocity[0] = velocity_temp[0];
                p1_new.velocity[1] = velocity_temp[1];
                p1_new.velocity[2] = velocity_temp[2];
                p1_new.position[0] = p1.position[0];

                //create the ionized ion (species3)
                cell_particles3->add_particles(1);
                id_new = cell_particles3->size() - 1;
                auto& p3_new = cell_particles3->data_[id_new];
                p3_new.velocity[0] = p2.velocity[0];
                p3_new.velocity[1] = p2.velocity[1];
                p3_new.velocity[2] = p2.velocity[2];
                p3_new.position[0] = p2.position[0];
                collide_times++;

                i_space_global = params->local_min[0] / params->cell_length[0] + i_space;
                diag1D->radiative_energy_collision[collision_number][i_space_global] += energy_threshold;
            }
        }
        cell_particles2->delete_useless_particles();

    }
}



